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Abstract 

We study the self-similar magnetohydrodynamics (MHD) of a quasi-spherical expanding void (viz. cavity or bubble) in the centre 
of a self-gravitating gas sphere with a general polytropic equation of state. We show various analytic asymptotic solutions near 
the void boundary in different parameter regimes and obtain the corresponding void solutions by extensive numerical explorations. 
We find novel void solutions of zero density on the void boundary. These new void solutions exist only in a general polytropic gas 
and feature shell-type density profiles. These void solutions, if not encountering the magnetosonic critical curve (MCC), generally 
approach the asymptotic expansion solution far from the central void with a velocity proportional to radial distance. We identify 
and examine free-expansion solutions, Einstein-de Sitter expansion solutions, and thermal-expansion solutions in three different 
parameter regimes. Under certain conditions, void solutions may cross the MCC either smoothly or by MHD shocks, and then 
merge into asymptotic solutions with finite velocity and density far from the centre. Our general polytropic MHD void solutions 
provide physical insight for void evolution, and may have astrophysical applications such as massive star collapses and explosions, 
shell-type supernova remnants and hot bubbles in the interstellar and intergalactic media, and planetary nebulae. 
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1. Introduction 

Supernova explosions, planetary nebulae and stellar winds 
from massive stars are believed to be the main sources of cre- 
ating voids (i.e. bubbles, cavities) in the interstellar medium 
(ISM) (e.g. Ferriere 1998, 2001 and extensive references 
therein). In the local ISM, the remnant of a typical isolated 
supernova grows for ~ 1.5 Myr and reaches a maximum radius 
of ~ 50 pc. The shell-type supernova remnants (SNRs) appear 
quasi-spherical and the masses within them have been swept up 
by ejecta from supernovae. An example of such structure is a 
void region towards the Lupus dark cloud complex, which was 
revealed through observations of 100 //m emissions (e.g. Gahm 
et al. 1990; Franco 2002) soft X-ray (e.g. Riegler et al. 1980), 
and 21 cm HI line (e.g. Colomb et al. 1984). Neutral hydro- 
gen (HI) voids have also been found in filled-centre SNRs (e.g. 
Wallace et al. 1994). 

Voids may also emerge long before we actually observe 
SNRs. According to the neutrino-driven mechanism of type- 
II and type-Ibc supernova explosions (e.g. Janka & Hille- 
brant 1989; Janka & Muller 1995, 1996), within the first sec- 
ond of a type II supernova, the intense neutrino flux generated 
by the central core bounce heats the surrounding stellar mass 
and pushes the stellar material outwards. Subsequently, a re- 
bound shock emerges and propagates outwards (e.g. Lou & 
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Wang 2006, 2007; Hu & Lou 2009). After several hundred 
milliseconds, the neutrino-sphere decouples from the gas, and 
may leave behind a cavity around the centre during a supernova. 
After this decoupling, the exploding star with a central cavity 
continues to expand and the central cavity eventually evolves 
into a hot bubble in the ISM. 

Massive stars can also create voids in the ISM through pho- 
toionization heating and stellar winds (e.g. Castor et al. 1975; 
Weaver et al. 1977; McKee et al. 1984). Likewise, fast winds 
from central compact hot white dwarfs may generate expanding 
cavities in planetary nebulae (e.g. Lou & Zhai 2009). HI voids 
and shells have been found and observed by radio observations 
around several Galactic Wolf-Rayet (WR) stars - massive stars 
undergoing significant mass losses (e.g. Cappa & Miemela 
1984; Cappa et al. 1986, 1988; Dubner et al. 1990; Niemela 
& Cappa 1991; Arnal & Mirabel 1991; Arnal 1992). On larger 
scales, observations show that voids also exist in neutral hydro- 
gen discs of spiral galaxies (e.g. Crosthwaite & Turner 2000). 
Recently, central cavities of ~ 200 kpc diameter and large-scale 
shock fronts have been revealed by Chandra X-ray observations 
in the galaxy cluster MS0735.6+7421 (e.g. McNamara et al. 
2005). 

The dynamic evolution of voids in the ISM still lacks a sys- 
tematic theoretical exploration. Dyson & Williams (1997) pro- 
vided qualitative description on the gas dynamic effects of mas- 
sive stars on the ISM. Chevalier (1997) studied the expansion 
of a photoionized stellar wind in late stages of stellar evolution 
(e.g. supernovae and planetary nebulae; see also Meyer 1997). 
From the centre, a stellar system consists of a hot bubble of 
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shocked fast wind, a region of shocked and photoionized wind, 
and an outer region of slow wind. Chevalier (1997) also em- 
ployed the isothermal self-similar transformation as Shu (1977) 
but without gravity to model the self-similar dynamic evolution 
of an outer slow wind. Physically, the fast hot wind bubble re- 
sembles the concept of a central void of this paper. Hu & Lou 
(2008a) presented self-similar void solutions to model "cham- 
pagne flows" of H II regions after the nascence of a massive pro- 
tostar in a conventional polytropic gas. Such voids embedded in 
nebulae can be created and sustained by fast stellar winds and 
photoionization heating. Here, we formulate an MHD problem 
with a general polytropic gas under self-gravity. 

The system of interest is a general polytropic magnetofluid 
with a quasi-spherical symmetry under self-gravity, thermal 
pressure gradient force and magnetic Lorentz force. We fur- 
ther find that MHD shocks are indispensable to establish sensi- 
ble global solutions, for example with the asymptotic velocity 
at large radii tending to zero. Magnetic field can be extremely 
important in many astrophysical processes on different scales 
and in particular, for star formation activities at various stages 
(e.g. Shu et al. 1987; Myers 1998). The Crab Nebula is ob- 
served to be supported by the magnetized pulsar wind (e.g. Lou 
1993; Wolf et al. 2003). When rotation is sufficiently slow in an 
astrophysical system, the overall geometry may remain quasi- 
spherical and the quasi-spherical random-field approximation 
(e.g. Zel'dovich & Novikov 1971) can be applicable. Chiueh & 
Chou (1994) discussed the gravitational collapse of an isother- 
mal magnetized gas cloud by including the magnetic pressure 
force from a randomly distributed magnetic field. Recently, we 
have provided detailed analyses by assuming a random trans- 
verse magnetic field with the consideration of both magnetic 
pressure and tension forces (Yu & Lou 2005; Yu et al. 2006; 
Lou & Wang 2007; Wang & Lou 2007, 2008). We presume 
that the non-spherical flows as a result of the magnetic tension 
force may be neglected as compared to the large-scale mean ra- 
dial bulk motion of gas. The key point is that the (self-)gravity 
is strong enough to hold on the entire gas mass and induce core 
collapse, or the driving force is strong enough to generate a 
quasi-spherical void. Therefore on large scales a completely 
random magnetic field contributes to the dynamics in the form 
of the average magnetic pressure gradient force and the average 
magnetic tension force in the radial direction. This approxima- 
tion was discussed in more details by Wang & Lou (2007). 

The general equation of state is p = K{r, t)p y where p is the 
thermal gas pressure, p is the mass density, y is the polytropic 
index and k is a proportional coefficient dependent on both ra- 
dius r and time t. For a global constant k, the equation of state 
is that of a conventional polytropic gas (e.g. Suto & Silk 1988; 
Lou & Gao 2006; Lou & Wang 2006, 2007; Hu & Lou 2008a). 
By setting y = 1 and k as a global constant, a conventional poly- 
tropic gas then becomes an isothermal gas. In case of y > 1.2, 
novel quasi-static asymptotic solutions for a polytropic gas ex- 
ist in approach to the system centre (see Lou & Wang 2006). 
As the specific enthalpy is p/(y - 1), we thus require y > 1 
to ensure a positive specific enthalpy. A general polytropic 
gas features the conservation of specific entropy along stream- 
lines. The conventional polytropic case is a only special case 



with constant specific entropy everywhere at all times. A gen- 
eral polytropic model with random magnetic field is the most 
general model of a polytropic magnetofluid of quasi-spherical 
symmetry under self-gravity (Wang & Lou 2008; Jiang & Lou 
2009). 

Our self-similar transformation employs a dimensionless in- 
dependent similarity variable x defined as a combination of ra- 
dius r and time t such that x = r/(k 1/2 t n ) where k is the so-called 
'sound parameter' to make x dimensionless and n is a key scal- 
ing index. Theoretically, our void solutions are those solutions 
whose enclosed mass is zero within a certain radius denoted as 
x*. This radius expands with time in a self-similar manner, i.e. 
r* = k l/2 x*t n . By this expression, the physical meaning of the 
self-similar scaling index parameter n is evident. The expan- 
sion speed of the void boundary is u* oc f (n_1) ; therefore for 
n > 1, a void expands faster and faster (i.e. acceleration), while 
for n < 1, a void expands slower and slower (i.e. deceleration); 
and for n = 1, a void expands at a constant speed. The void 
boundary may also be regarded as an idealization of a contact 
discontinuity between a faster wind and a slower winds (e.g. 
Chevalier 1997; Lou & Zhai 2009). 

The case of y = 4/3 corresponds to a relativistically hot gas 
that deserves a special attention. Homologous core collapse 
for a relativistically hot gas was studied by Goldreich & Weber 
(1980) and the behaviour of such system has been treated by 
Yahil (1983) as a limit of y -> (4/3) + . Recently, Lou & Cao 
(2008) presented an illustrative example of void in such a sys- 
tem. In this paper, we study voids for a magnetized Newtonian 
gas (i.e. y + 4/3 in general), as well as voids in a relativistically 
hot fluid (i.e. y = 4/3), and we offer several concrete examples. 

This paper is structured as follows. Section 1 is an introduc- 
tion to provide background information and the motivation of 
this investigation; Section 2 describes first the formulation of a 
general polytropic magnetofluid under quasi-spherical symme- 
try, and secondly analytic asymptotic solution behaviours near 
the void boundary in various parameter regimes, and thirdly 
various asymptotic MHD solutions that are useful in construct- 
ing global semi-complete solutions (i.e. solutions that are valid 
in the range < x < +°o); Section 3 describes and discusses 
properties of void solutions with different parameters, in con- 
texts of hydrodynamics and MHD, and presents a few exam- 
ples; Section 4 gives examples of astrophysical applications of 
such void solutions. Finally, Section 5 contains conclusions 
and discussion. Mathematical derivations are included in an 
appendix. 

2. Self-Similar MHD with Quasi-Spherical Symmetry 

2.1. Formulation of a Nonlinear MHD Problem 

The MHD evolution of a general polytropic gas of quasi- 
spherical symmetry and under self-gravity can be described by 
a set of nonlinear MHD partial differential equations (PDEs) in 
spherical polar coordinates (r, 6, </>), namely 
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where p(r, f) is the mass density, u(r, t) is the bulk radial gas 
flow speed, M(r, t) is the enclosed mass within radius r at time 
t, p is the thermal gas pressure, G = 6.67 x 10~ 8 dyne cm 2 g~ 2 is 
the gravitational constant, and < B 2 > is the ensemble average 
of a random transverse magnetic field squared (i.e. proportional 
to the magnetic energy density). Equations (HJ and d2j represent 
mass conservation, leading to dM /dt+udM /dr = 0. We assume 
a random magnetic field mainly in transverse directions, and the 
magnetic force perpendicular to magnetic field lines directs to 
the radial direction and appears in the radial momentum equa- 
tion ® as the magnetic pressure and tension terms. Together 
with magnetic induction equation 
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with y being the polytropic index, we complete the model for- 
mulation for a general polytropic MHD with a quasi-spherical 
symmetry (Wang & Lou 2007, 2008). 

In this paper, we consider self-similar solutions which form 
an important subclass of nonlinear MHD PDEs. In order to re- 
duce these nonlinear MHD PDEs to self-similar ordinary dif- 
ferential equations (ODEs), we introduce the following self- 
similar transformation as Wang & Lou (2008), namely 
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where v(x), a(x), B(x), m(x), w(x) are dimensionless reduced 
variables of x only. We refer to v(x), a(x), B(x), m(x) and w(x) 
as the reduced radial speed, mass density, thermal pressure, en- 
closed mass and magnetic energy density, respectively. Self- 
similar transformation (O is identical with that of Wang & Lou 
(2007, 2008). By substituting self-similar transformation © 
into equations (T])-©, we obtain several valuable integrals 



w - ha 2 x 2 , 
B = Ca 7 m q , 
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m — ax l (nx — v) . (9) 

Equation © corresponds to the frozen-in condition for mag- 
netic field in the ideal MHD approximation, where h =< B 2 t > 
l(\6n 2 Gp 2 r 2 ) is a dimensionless magnetic parameter represent- 
ing the average strength of a random transverse magnetic field. 

Equation ([8]) is the self-similar form of the specific entropy 
conservation along streamlines, where the exponent parameter 



q = 2(n + y - 2)/(3n - 2) and thus y — 2 - n + (3n - 2)q/2, 
and C is an arbitrary coefficient from integration. For q — 0, 
the flow system involves a conventional polytropic gas with a 
constant specific entropy everywhere in space and at all times; 
for q > 0, the specific entropy increases from inside (smaller 
x) to outside (larger x); and q = 2/3 leads to y — 4/3 for a 
relativistically hot gas (e.g. a photon gas, a neutrino gas or an 
extremely high temperature electron gas) with an arbitrary c. 
Actually we may set C = 1 in all cases with q + 2/3 without 
loss of generality, because an adjustment of sound parameter k 
in self-similar transformation (O to C l ^ l ~ iql2) k would make C 
disappear. 

Equation (0 requires both 3« - 2 > and nx - v > for a 
positive enclosed mass M(r, t). When nx— v = at a certain x*, 
the enclosed mass vanishes by equation (0; this is referred to as 
a void with x* being the void boundary in a self-similar expan- 
sion. Accordingly, the reduced radial velocity on the boundary 
is given by v* = nx*. The condition nx - v = on the void 
boundary implies that the self-similar expansion speed of the 
void boundary dr* /dt is equal to the radial flow velocity on the 
void boundary u(r* ', t). This is regarded as the physical condi- 
tion for a contact discontinuity between the outer slower stellar 
wind and the inner faster wind driving a hot bubble by Chevalier 
(1997) for q -Q,n - 1, y = 1 and without gravity. Lou & Zhai 
(2009) considered a gas dynamic model for planetary nebulae 
with contact discontinuities for an isothermal self-gravitating 
gas. From now on, we denote variables on the void boundary 
by a superscript asterisk *. 

Combining all equations above, we obtain coupled nonlinear 
MHD ODEs for the two first derivatives a' and v' in the follow- 
ing compact forms of 



X{x, a, v)a' — J[{x, a, v) , 



X(x,a,v)v' =<V(x,a,v) , (10) 



where the three functional coefficients X, 31 and <V are explic- 
itly defined by 
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The formulation above is largely the same as Wang & Lou 
(2008), except for an additional free parameter C in cases with 
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q = 2/3 (i.e. y = 4/3). Wang & Lou (2008) also provides 
procedures to determine magnetosonic critical curve (MCC), 
eigensolutions across the MCC and MHD shock jump condi- 
tions across the magnetosonic singular surface (see Appendix 
EJ. 

The MCC for q — 2/3 and thus y — 4/3 is special. Extensive 
numerical explorations suggest that the MCC bears the simple 
form of v = Tjx, a = constant, where 77 is a constant coefficient 
dependent on parameters n, h and proportional factor C. Sub- 
stituting this form into equation (TTTT) and the MCC conditions 
X — and become respectively 



a — 
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h + (4/3)C(n - j]) 2 ' 3 ' 



Cain - 7?r I/3 (2 - 3n) + 2(1 - t])(n - 77) 
= (n - 77)77 + + 2ha . 
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Substituting relation (fT2l into equation ([TBI to eliminate a, we 
immediately derive an expression for the constant coefficient 77 
in terms of 77, h and C. Once 77 is known in relation (fT2l . we can 
compute a value accordingly. With /1 = 0, relations (fT2l and 
( fT3l give the same solution as discussed in Lou & Cao (2008) 
for a nonmagnetized relativistically hot gas. Here, we extend 
the special MCC with y = 4/3 to a magnetofluid embedded 
with a random transverse magnetic field. 

2.2. Behaviours of a Polytropic Void Boundary 

We now analyze asymptotic behaviours around the void 
boundary x*, and these boundary conditions will be used to con- 
struct various solutions in numerical integrations starting from 
the void boundary. 

The gas pressure should be continuous across the void 
boundary, otherwise a shrinkage of the void boundary with dif- 
fusions would be expected. According to equation © and for 
a + 0, inequality q < on nx* - v — leads to a diverging re- 
duced pressure B. For q > 0, it follows automatically that B = 
at the void boundary. Therefore to ensure B* = B(x*) = 0, we 
require a* = a(x*) = in cases of q < 0, or a* > in cases of 
q > 0. It is favorable to further require the continuity of mass 
density, as a discontinuous density would lead to a local diffu- 
sion, in addition to a global self-similar evolution. Therefore, 
a solution with a = at nx* - v = is regarded as a physi- 
cally sensible one, otherwise the self-similar solution should be 
seen as an asymptotic solution valid in the region sufficiently 
far from the void boundary. 

2.2.7. Hydrodynamic and MHD Cases with a* — 
on the Expanding Void Boundary 
The void boundary obeying a — and nx - v — may pos- 
sibly become a critical curve with the three functional coeffi- 
cients on both sides of equations ( TTOb being zero. According 
to equation ( fTTT i. the possible non-zero terms of the three func- 
tional coefficients X, and *V approaching the void boundary 



77.x - v = and a — are the thermal pressure gradient force 
terms, namely 

X~Cyx 2q a y+q - x (nx-v) q , 
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<V ~ C[2(l - n)y + (2 - 3n)q]x 2q a y+q - 1 (nx - v) q 
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where we have used the relation y — 2 - n + (3n - 2)q/2 . 

According to expressions (TT~4-b . the parameter regime of q > 
, y + q > 1, except for the special isothermal case (q = 0, y + 
q — 1), ensures the vanishing of X and *V on the void boundary; 
when a ~ (nx - v) near the void boundary as shown by our 
analysis presently, then J[ also vanishes at the void boundary, 
and the void boundary indeed becomes a MCC. 

By a local first-order Taylor series expansion, we obtain from 
equations (TTOb and (fTTT i two pairs of eigensolutions for the first 
derivatives of v(x) and a(x) across the void boundary as a criti- 
cal curve. The one that ensures positive enclosed mass is 
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It can be shown that for 77 + 1, all points on this critical curve 
are saddle points (e.g. Jordan & Smith 1977). It is known 
that around a saddle singular point, only solutions along the 
direction defined by eigensolutions are allowed. Therefore for 
a* — at the void boundary, the self-similar solutions would 
follow the behaviour described by expression (fTsl l. The pres- 
ence of magnetic field is crucial here. In a purely hydrodynamic 
case of h = or a conventional polytropic case of q — 0, the re- 
duced density gradient a' diverges at the void boundary x — x*. 
It would require additional considerations for a void boundary 
on which the mass density vanishes but its first derivative di- 
verges. To better describe a sudden change of mass density on 
the void boundary, we would then set a non-zero mass density 
a* + at the void boundary (see analyses in the following sec- 
tions). The valid regime of parameters in which solution (TT~5T > 
stands is therefore g>0, y + q > 1 and h > 0. As we should 
further require y > 1 for a positive enthalpy, the valid regime of 
parameters is simply q > and h > 0. 

If the three functional coefficients X, J[ and 'V do not vanish 
at the void boundary, the leading terms of the first derivatives 
of v and a at the void boundary are then 



v' |. v .= 2(1 -«) + (2 -3n)q/y : 
(2 - 3n)q a 
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Note the coefficient C disappears in these two expressions. The 
second equality for a' in equation ( fTTh involves equation ( IT6b . 
The asymptotic solution of a near the void boundary is 

a = K(x - x*)- qKy+q) + • • • , (18) 
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Figure 1 : Regimes of parameters q and y for which asymptotic solutions near 
the void boundary are applicable. The horizontal axis is the parameter q and the 
vertical axis is the polytropic index y. The solid line is for y+q = 0. The vertical 
axis of q = corresponds to the conventional polytropic case (n + y = 2) and 
the point q = Q , y = \ (n = \) corresponds to the isothermal case. The light 
shaded zone in the second quadrant is the regime to apply asymptotic solution 
i 1 6i - i 1 8i and the heavy shaded zone in the first quadrant is the regime to 
apply asymptotic solution \\5\ . We only draw the y > part out of physical 
consideration; in fact, we should always require y > 1 to ensure a positive 
specific enthalpy. 



where K is an arbitrary integration constant referred to as the 
void density parameter. To ensure solution < TT~8T > going to zero 
in the limit of x — > x*, we require q/{y + q) < 0. We now 
verify that for such solutions, X, 3i and *V actually diverge at 
the void boundary. With solutions (fT6]>-(fT8b. we have X ~ 
(x- X y l(y+c '\<V ~ (x- X y l(y+c '\jl ~ (x-x*y l . Therefore, the 
condition q/(y + q) < in the meantime ensures the validity of 
the asymptotic solution. We should require y > 1 as a physical 
requirement and thus inequality q/(y + q) < is equivalent to 
q < and y + q > 0. 

A positive enclosed mass requires v' | A .< n, which provides 
a lower limit for the index q, namely 



q > q n 



4 + 3n 2 - Sn _ 2(n - 2) 
3n(3n/2- 1) ~ ~ln 



(19) 



where n — 2 < 0. With algebraic manipulations, it is proven that 
inequality ( fT9l together with y > 1 is equivalent to inequality 
y + q > 0. Thus within the parameter regime of q < and 
y + q > where solutions (fTol l - ( fToT l stand, inequality v' \x*<n 
is automatically satisfied. 

In summary, in the vicinity of the void boundary nx - v — 
with a — 0, we find two possible types of asymptotic solutions. 
The parameter regimes in which these solutions are applicable 
are shown in Figure Q] accordingly. With q > (the heavy- 
shaded zone in the first quadrant of Figure Q3 and h > (mag- 
netized), the asymptotic solution takes form ( fT5l >. referred to as 
LH1 solutions. With q < 0, y + q > (the light-shaded zone 
in the second quadrant of Figure [TJ, the asymptotic solution 
takes form ( fTol l - (fTol l. referred to as LH2 solutions. For other 
regimes of parameters, no void solutions with boundary condi- 
tion nx-v - and a — are found. The isothermal case and the 
conventional polytropic case do not satisfy either of these two 
requirements, so they need additional considerations for asymp- 
totic solutions near the void boundary with a* + (e.g. Hu & 



Lou 2008a; Lou & Zhai 2009). Hence, the novel LH1 and LH2 
void solutions only exist in general polytropic MHD cases with 
q > and q < 0, respectively. LH1 void solutions require the 
presence of magnetic field, while LH2 solutions remain valid in 
the purely hydrodynamic case as well. 

2.2.2. Hydrodynamics of a* > at the Void Boundary 

The situation with a* + at the void boundary is intrinsically 
different. In such cases, we require q > to ensure the pressure 
going to zero at the void boundary. One can then numerically 
integrate coupled nonlinear ODEs (TTOt and (fTTb directly from 
the void boundary to construct solutions. The possible singular- 
ity at the void boundary, and the corresponding leading terms 
for X, 3i and *V approaching the void boundary nx - v — 
depend largely on parameter q. We find that q = for a con- 
ventional polytropic gas is a special case giving a different ex- 
pression of the asymptotic behaviour at the void boundary. We 
examine below two situations of q = and of q > separately. 
Case of q = 

In the q — case, the formulation is simplified considerably 
and we return to the conventional polytropic case of n + y = 
2. For non-magnetic cases of h — 0, the asymptotic solution 
approaching the void boundary is 



nx* + 2(1 — n)(x — x*) + ■ 



a — a + 



n(l - n) 



a x (x — x ) + ■ ■ 



(20) 



(21) 



where x* denotes the void boundary in a self-similar expansion 
and a* denotes the reduced mass density on the expanding void 
boundary. This solution is the same as the void solution in Hu 
& Lou (2008a). Expression (|2TT> hints that, if a* = 0, then 
the solution becomes a(x) = everywhere for all x, and this 
is indeed so. Again, a* — is not allowed for a conventional 
polytropic gas. In this case, no apparent singularity is found 
near the void boundary, and both terms (nx - v) and (a - a*) 
scale as (x - x*), which we denote as type-N (Normal) void 
behaviour. Note that LH1 void solutions represent also a type- 
N void behaviour. 

Series expansion solutions (l20l i and d2Tb become insufficient 
for n — 1. In this isothermal case of q = and n — 1 (y — 1), 
we can obtain asymptotic solutions near the void boundary to a 
higher order. The leading terms of v' as x — > x* yield 



v' \ x *= — (x-v) 



We then obtain the leading terms of v and a as x — > x* 



* 1 » 2 

V — X H (x — X ) + ■ 

X* 



x*f + ■ 



(22) 



(23) 



(24) 



No singularity appears in the isothermal case (as the Type-N 
behaviour) but the term {nx - v) has the leading order of magni- 
tude (x - x*) 2 , which we refer to as the type-N2 void boundary 
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(see Lou & Zhai 2009 for isothermal voids in self-similar ex- 
pansion). 

Cases of q > 

In such cases, the leading terms of functional coefficients X, 
J{ and 'V are also the thermal pressure terms as in equation (fT4l i. 
and then v' | x « and a' | x . are the same as expressions (TTSl l and 
( fTTI i. Therefore, the asymptotic form of a approaching the void 
boundary is the same as equation ( [TBI , viz. a = K(x-x*)~ q ^ r+q) 
where K is the void density parameter. The first derivative of v 
in these cases tends to a certain value on the void boundary and 
the term nx - v has the leading order of magnitude x - x*. 

As q > in this case, a diverges on the void boundary. 
Here a sharp discontinuity in mass density exists across the void 
boundary. We verify in turn that the divergence of a does not 
affect the leading terms of functional coefficients X, J[ and *V 
and the validity of expressions ( fT6b and dTTb . Despite this diver- 
gence, m(x) remains continuous at x = x* (equation[9]) and the 
asymptotic solution with such singularity is physically allowed. 
We refer to such asymptotic solution on the void boundary as 
the type-D (diffusion) void behaviour. 

We may regard the void boundary as a translation of centre 
along the streamline nx - v — 0. Previous asymptotic solutions 
at x — give either a zero a or a divergent a obeying power 
law (e.g. a oc x~ 3 ^ 2 for free-fall solutions, see Lou & Wang 
2007). Here on a void boundary, the power law index of the 
asymptotic a depends on parameter q. In terms of physics, a 
local diffusion process may smooth out this singularity, bear- 
ing in mind that self-similar behaviours will be modified by lo- 
cal diffusions near the void boundary in self-similar expansion. 
Relevant comments on this may be found in Lou & Cao (2008) 
and Lou & Zhai (2009). 

With q = in expressions ( TToT ) and ( fTTT i. we have the same 
v' \x> as equation ( f20b and a* = K. In fact, the q — case 
(Type-N; a tends to a positive constant) is a transitional case 
between q < (LH2, a tends to zero at the void boundary) and 
q > (Type-D; a diverges at the void boundary). So far we 
have provided sensible void solutions for purely hydrodynamic 
cases with different q values. For all such solutions, the ther- 
mal pressure force becomes dominant near the void boundary. 
Without magnetic field, the thermal pressure is the key factor in 
determining the dynamics near the void boundary as it should 
be. 



expressions (f20b and (f2TT) . 



v — nx + - — — : — — — (jc— x ) + ■ 



a - a - 



+ ha*(x*) 2 

n(n- IK + 2h(a*f 
y(a* )('-") + ha*(x*) 2 



x*(x - x*) + ■ ■ ■ 



(25) 
(26) 



Setting h — in solutions d25l ) and d26l l. we retrieve solutions 
(l20b and ( |2"TT i as a necessary check. This MHD solution mani- 
fests a type-N void behaviour and the type-N2 void behaviour 
of the isothermal case disappears. 
Cases of < q < 1 

In such cases, the presence of magnetic field becomes the 
leading term approaching the void boundary, and the first 
derivatives of v and a are respectively 



2(1 -n) 

a* 



(27) 



a' \x'= C 



?(2-3n) ,,.i_„ 



h 



(or') 1 



" +3 ^V) 2<? - 2 (nx - v)i- 1 . (28) 



Equation ( I28l l gives the leading term of the asymptotic solution 
of a near the void boundary as 



a — a + C 



— (2 -^(aT(**) 2 - 2 



(29) 



x(a*n + 2n-2) q -\x- x*y 



To ensure the validity of this solution, we should require v' |.^< 
n for a positive mass. For n > 1, this condition is satisfied 
automatically, while for n < 1, the condition implies a*n + 2n- 
2 > or equivalently n > 2/(2 + a*). No apparent singularity 
exists in this solution and the term (a - a*) scales as (x - x*) q . 
We refer to this asymptotic solution as type-Nq void behaviour. 
Cases with q = 1 

The asymptotic behaviour of v near the void boundary is the 
same as equation ( l27l i. while the first derivative of a becomes 



(a*) 1+n/2 x*(2 - 3n) - (n - l)n - 2ha* 
hx* 



(30) 



2.2.3. MHD Cases with a* + at the Void Boundary 

In such cases, we first require q > to ensure the pressure 
approaching zero at the void boundary. With h > in equations 
dTOb and ( fTTl ). the leading terms in the vicinity of a void bound- 
ary is different from the hydrodynamic case. In the presence 
of magnetic field, the magnetic force becomes dominant at the 
void boundary and diffusion behaviours (Type-D) of the void 
boundary do not appear. The void boundary generally shows 
no singularity in these MHD cases. There are four distinct situ- 
ations described below in different regimes of q parameter. 
Case q = for a conventional polytropic gas 
In such cases of conventional polytropic MHD, additional 
terms associated with the magnetic force appear in the asymp- 
totic solution, having a similar form parallel to hydrodynamic 



giving a type-N behaviour near the void boundary. For h — » 0, 
the first derivative a' would diverge. 
Cases with q > 1 

The asymptotic behaviour of v near the void boundary re- 
mains the same as equation d27l ). while the first derivative of a 
becomes 

(n-l)n + 2ha* 



a x«= -- 



hx* 



(31) 



again giving a type-N behaviour near the void boundary. 

In different parameter regimes of q, n (or y) and h, we have 
obtained different types of asymptotic behaviours near the void 
boundary and different expressions of asymptotic solutions. We 
summarize these results in Table [T]for reference. 
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Table 1 : Summary of asymptotic solution behaviours near the self-similar void boundary nx — v = and a # 0. In each case, we show the behaviour type and then 
equation numbers of the corresponding asymptotic solution. The properties of different types of solutions on the void boundary are described below the table. 



<7 


h — 

n — \j 


n j- u 


q = 


Type-N, d20l|2B; Type-N2 <|23]|2l 


Type-N, $Z5\UB 


0<q<l 


Type-D, dHHUi 


Type-Nq, (|27l|28]) 


q = l 


Type-D, (HBE) 


Type-N, (HZl [30j» 


q > 1 


Type-D, ([161 [17} 


Type-N, (|27l [3TJ 



• Type-N: a tends to a nonzero finite value; and (a - a") and (v - v*) scale as (x - x*). 

• Type-N2: a tends to a nonzero finite value; and {a - a*) and (v - v*) scale as (x - x*) 2 . 

• Type-Nq: a tends to a nonzero finite value; (a - a*) scales as (x - x*) q ; and (v - v*) scales (x - x*). 

• Type-D: a diverges and scales as (x - x*Y q ^ y+q) \ and (v - v*) scales (x - x*). 



2.3. Asymptotic Self-Similar Solutions at Large x 

Prior studies have revealed various asymptotic self-similar 
solutions of a quasi-spherical magnetofluid under self-gravity. 
At small x, we have derived quasi-magnetostatic asymptotic 
solutions (Lou & Wang 2006, 2007; Wang & Lou 2007 with 
q = 0; Wang & Lou 2008 with q + 0), central MHD free-fall 
solutions (Shu 1977 for an isothermal gas; Suto & Silk 1988 for 
a conventional poly tropic gas; Wang & Lou 2008 for a general 
polytropic gas), strong-field asymptotic MHD solutions (Yu & 
Lou 2005 for an isothermal gas; Lou & Wang 2007 for a con- 
ventional polytropic gas; Wang & Lou 2008 for a general poly- 
tropic gas). At large x, we have asymptotic MHD solutions 
described below. 

2.3.1. Asymptotic MHD Solutions of Finite Density 
and Velocity in the Regime of Large x 
In this case, the gravitational force, the magnetic force (i.e. 
the magnetic pressure and tension forces together), and the ther- 
mal pressure force are in the same order of magnitude at large 
x. The asymptotic solutions at large x are given by Wang & Lou 
(2008), namely 



a = Ax 



■2/n 



v = Bx 



1-1/n + | 



2h(n- 1) 



_(3n — 2) n 
+2(2 - n)n q - 1 A 1 -" +3 "i /2 \x 1 - z/n + • • 



(32) 



where A and B are two constants of integration, referred to as 
the mass and velocity parameters respectively. To ensure the 
validity of solution d32b . we require 2/3 < n < 2 (note that 
inequality n > 2/3 is directly related to self-similar transforma- 
tion © and a positive enclosed mass). In case of 2/3 < n < 1, 
the mass and velocity parameters A and B are fairly arbitrary. 
In case of 1 < n < 2, velocity parameter B should vanish to 
ensure that v tends to zero at large x. This valid range of scal- 
ing parameter n corresponds top cc r~ 3 to p oc r _1 . For the 
dynamic evolution of protostellar cores in star-forming clouds, 
power-law mass density profiles should fall within this range. 



This appears to be consistent with observational inferences so 
far (e.g. Osorio, Lizano & D'Alessio 1999; Franco et al. 2000; 
McKee & Tan 2002). 

Furthermore by setting v = in MHD ODEs ([TOj and (flTl 
we readily obtain an exact global solution in a magnetostatic 
equilibrium, namely 



-2/n 



■ AqX 



where the proportional coefficient Aq is given by 



n z - 2(2 - «)(3« - 2)h 



2(2 - n)(3n - 2) 



-l/(«-3«ij/2) 



(33) 



(34) 



This describes a more general magnetostatic singular polytropic 
sphere (SPS) with a substantial generalization of q ^0; the 
case of q = or n + y — 2 is included here and corresponds 
to a conventional polytropic gas of constant specific entropy 
everywhere at all times (Lou & Wang 2006, 2007; Wang & 
Lou 2007, 2008). 

2.3.2. Asymptotic MHD Thermal Expansion Solutions 

At large x, the pressure may become dominant in certain sit- 
uations (Wang & Lou 2008). Thus we may drop the magnetic 
and gravity force terms in ODEs ( TTOb and (fTTT i. By assuming 
v ~ cx + b and a ~ Ex p with c, E and P being three constant 
coefficients, ODEs ([Toll and (Qj} then lead to 



P = - 



Oq-2) 



(l-n + 3nq/2) ' 



?l-n+3nq/2 



(n-c) q (2 + P) = c(l -c) 



P = 



(3c -2) 
(n-c) 



(35) 



where the three constant coefficients (c, E, P) can be deter- 
mined by equation d35t . This solution for MHD thermal ex- 
pansion is valid for q > 2/3 as we need a power-law exponent 
P < for a converging a(x) at large x. We note that for a cer- 
tain system whose parameters are predefined, only one thermal 
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expansion solutions at large x is allowed, except for a free pa- 
rameter b. Actually, a translation on v will not alter the structure 
of the solutions. The radial bulk flow speed at large x is 



(36) 



At a certain time t, the radial flow speed is simply proportional 
to r. There is a qualitatively similar flow speed profile in the 
special case of y — 4/3 for a relativistically hot gas (Goldreich 
& Weber 1980; Lou & Cao 2008; Cao & Lou 2009). 



2.3.3. MHD Free-Expansion Solution 

In cases of q < 2/3, numerical exploration suggests an ex- 
pansion solution in the asymptotic form of 



2 

—x + 1 

3 



+ O0 , 



(37) 



where is a constant value of a at large x. With the radial 
velocity proportional to the radius in asymptotic form (f3Tb and 
q < 2/3, the pressure gradient terms in nonlinear MHD ODEs 
( [Tol l and (fTTT > can be dropped and the leading terms of the three 
coefficients X, and *V as x — > +oo are 



X ~ hax - (nx - v) 



x — V 

J{ ~ 2 a(nx — v) — a 

x 



(nx - v) 
(n — l)v + — —a + 2nax 



'V ~ 2- -/zoo: 2 — (nx — v) 

x 



(3n - 2) 



(nx — v) 
(n — l)v + — —a + 2hax 



(3n - 2) 



(38) 



respectively. To obtain asymptotic solution (l37l >. we require 
J{ — > and 'V/X — * 2/3. The constant then obeys the 
following relation 



(1 +6h)a 2 -2al3 = 0. 



Equation (1391 has only one non-trivial solution, namely 



3(1 + 6/7) 



(39) 



(40) 



We substitute this a x into condition 'V/X — > 2/3 and find that 
this condition is satisfied. The other solution — is indeed 
trivial and does not satisfy condition 'V/X — * 2/3. In sum- 
mary, we verify the existence of asymptotic expansion solution 
( 137] ) in the regime q < 2/3 and the constant value q- m is given 
by equation (l40b . For such an expansion solution, the pressure 
gradient is negligible, we thus refer to such expansion solution 
as the 'free-expansion' solution. The 'free-expansion' solution 
is the counterpart of thermal expansion solution for q < 2/3, as 
shown in this paper. The constant a x depends not on parameter 
n, but on magnetic parameter h. With a larger h (i.e. a stronger 
magnetic field), the constant asymptotic density is lower. 



2.3.4. The MHD Einstein-de Sitter Solution 

There exists a special exact semi-complete global solution 
referred to as the MHD Einstein-de Sitter solution, having the 
form of v = 2x/3 and a = constant for all x. Wang & Lou 
(2007) described this solution for the case of a conventional 
polytropic magnetofluid (i.e. q = 0). The form of this MHD 
Einstein-de Sitter solution is 



v = —x . 
3 



3(1 + 6/7) 



/ 2\ 2jc 3 

r " 3J3717 



3(1 +6/1) 



q = Q. 



(41) 



where n > 2/3. Compared with free-expansion solution dTTI i 
and ffco value ( l40b , we find that the 'free-expansion' solution 
naturally becomes the MHD Einstein-de Sitter solution with 
q — 0. We extend the consideration to the general polytropic 
form adopted in this investigation. By setting v = 2x/3, a = 
const in nonlinear MHD ODEs dTOb and (fTTT l. it is clear that 
only q = 2/3 (i.e. y = 4/3 with an allowed range of n), in addi- 
tion to the case of q = 0, will make the solution valid for all x. 
We have the novel MHD Einstein-de Sitter solution as 



v=-x. 



2 

oc — — 



1 



(n 



3 6A+1+ 6C(n 

2/3)x 3 



3 6/1 + 1+ 6C(n - 2/3) 2 ' 3 



■ 2/3) 2 / 3 ' 



q= 3 



(42) 



where n > 2/3. Comparing with equation (80) of Wang & Lou 

(2007) , we have a more general form of Einstein-de Sitter so- 
lutions for a relativistically hot gas, for which the gravity and 
magnetic forces cannot be neglected with respect to the pres- 
sure force. Lou & Cao (2008) studied a similar relativistically 
hot gas of spherical symmetry with another self-similar trans- 
formation and derive another form of Einstein-de Sitter solu- 
tion (36) of Lou & Cao (2008). With the freedom to choose 
C parameter, which is linked with parameter Co in Lou & Cao 

(2008) , these two forms of Einstein-de Sitter solution are equiv- 
alent for h = 0. It is interesting to observe that such globally 
exact Einstein-de Sitter solution can only exist either in a con- 
ventional polytropic gas or in a relativistically hot gas. 

We have explored above asymptotic expansion solutions in 
different q regimes and derived three kinds of expansion solu- 
tions. With q < 2/3, the thermal pressure can be neglected and 
the free expansion solution stands; with q > 2/3, the gravita- 
tional and magnetic forces can be neglected and the thermal- 
expansion solution stands; with q = 2/3, all the three forces 
are comparable and the Einstein-de Sitter solution is an exact 
global solution. Here we see that q = 2/3 separates different 
situations of self-similar expansion behaviour. All these expan- 
sion solutions have velocity proportional to the radius. For the 
free-expansion and Einstein-de Sitter solutions, v = 2x/3 and a 
is equal to a certain constant, while for the thermal-expansion 
solution, v - cx and a converges to zero at large x. We will 
see presently that in general MHD void solutions merge into 
one kind of expansion solutions (determined by q parameter) 
far from the flow centre, if the MCC is not encountered. 
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Figure 2: Global LH2 void solutions with shocks for parameters (n = 0.75, q = 
-0.5, y = 1.1875, h = 0) and the void boundary x* = 1. The upper panel shows 
the reduced density a(x) and the lower panel shows the reduced radial velocity 
v(x). In the upper panel, a logarithmic scale is used for a(x). In both panels, the 
dotted curve is the sonic critical curve, and in the lower panel the dashed curve 
is the void boundary line nx — v = 0. The void solution with K = 1 connects 
two solutions 1 and 2 whose parameters are: A = 7.3025, B = 3.7802, x sd = 
3, a si = 1.6366, v sd = 1.9702, x su = 3.6205, a sa = 0.2519, v su = 0.1559 
(solution 1), and A = 8.7808, B = 6.3723, x si = 3.2, a sd = 1.1040, v sd = 
1.9755, jc su = 3.2529, a su = 0.3781, v su = 1.1854 (solution2), respectively. 
The void solution with K = 5 connects two solutions 3 and 4 whose pa- 
rameters are: A = 24.493, B = 4.854, x si = 2.27, a sd = 2.8499, v sd = 
1.4515, x su = 4.5410, a su = 0.2587, v su = -2.1237 (solution 3), and 
A = 7.6878, B = 5.4448, x sd = 2.37, a sd = 2.0616, v sd = 1.4308, x w = 
2.5603, a su = 0.3928, v su = -0.0456 (solution4), respectively. 

3. Hydrodynamic and MHD Void Solutions with a* = 

3.1. Hydrodynamic Cases 

In purely hydrodynamic cases with h — 0, LH1 void solu- 
tions do not exist and we consider only LH2 void solutions. 
The parameter regime for LH2 void solutions is q < and 
y + q > 0. To illustrate LH2 void solutions by examples, 
we choose a set of parameters as (« = 0.75, q - -0.5, y = 
1.1875, h = 0) and construct LH2 void solutions with assigned 
values of void boundary x* and density parameter K of asymp- 
totic form dT6l> — (TT~8b. We choose a downstream shock position 
x sl j and insert a hydrodynamic shock there to match inner void 
solution with outer asymptotic envelope solution d32l of finite 
velocity and density at large x. Note that with n < 1 the veloc- 
ity actually tends to zero at large radii. Several such global LH2 
void solutions with shocks are shown in Figure [2] We have also 
performed numerical explorations with different parameter sets 
and the results are qualitatively similar. 

With a* — at the void boundary, the density first increases 
and then decreases as x increases; and the radial velocity in- 
creases as x increases. The density profiles of LH2 void so- 
lutions (see the upper panel of Figure |2]i indicate a prominent 



shell-type morphology surrounding a central cavity in expan- 
sion. The peak density of the solution and the width of the 
shell is modulated primarily by K parameter, which varies for 
different astrophysical gas flow systems. With different values 
of x s d, different dynamic behaviours of the corresponding up- 
stream side can be obtained (e.g. the lower panel of Figure |2). 
In the vicinity of the upstream shock front, the fluid can be ei- 
ther an inflow (solutions 1 and 2) or an outflow (solutions 3 
and 4). With adopted parameters, the fluid always merges into 
an asymoptotic outflow (parameter B > 0) far from the system 
centre, whereas it possible to construct solutions with B < for 
LH2 voids (see MHD examples below and examples in Figure 
[6] for a conventional polytropic gas). Numerical explorations 
indicate that it is generally not possible to make LH2 void solu- 
tions to cross the sonic critical curve smoothly. In other words, 
the inclusion of hydrodynamic shocks is necessary in order to 
construct sensible semi-complete global void solutions. 

The possibility of asymptotic inflows at large x associated 
with central void solutions deserves special attention. With- 
out shocks, void solutions generally merge into asymptotic ex- 
pansion solutions with flow velocities remaining positive. This 
establishes the physical links between central expanding voids 
and asymptotic outflows. In the presence of shocks, the up- 
stream may have inflows, either near the shock front or suffi- 
ciently far from the system centre. This means that initially 
and at large radii a cloud system may involve inflow or con- 
traction under the self-gravity. When the central engine forms 
an expanding void, a resulting shock may face the falling gas 
and expand outwards. For example, in Hu & Lou (2008a), the 
possibility of asymptotic inflows at large x is interpreted as a 
special scenario for "champagne flows" in H II regions. 

We emphasize that LH2 solutions here are the only void so- 
lutions in non-magnetized cases with a* = 0. The shell-type 
appearance is a general feature for LH2 void solutions. Such 
solutions are applicable to shell-type morphologies, widely ob- 
served in various astrophysical gas systems, such as supernova 
remnants and hot bubbles (e.g. Ferrierre 1998, 2001), H II re- 
gions (e.g. Hu & Lou 2008a) and even cavities in galaxy clus- 
ters (e.g. McNamara et al. 2005). The increasing velocity with 
radius of such solutions suggests a wind nature: the fast wind 
from the central cavity decelerates in the shell and the mass is 
accumulated in the shell. This is consistent with the picture of 
champagne flows of H II regions (e.g. Hu & Lou 2008a) and 
supernova remnants. 

3.2. MHD Cases 

With a random magnetic field, both LH1 and LH2 void so- 
lutions exist. As counterparts to hydrodynamic cases, we con- 
sider LH2 void solutions with the same parameters adopted in 
the previous section, except for the magnetic parameter being 
h - 0.3. Such global MHD LH2 void solutions with shocks are 
shown in Figure [3] 

The appearance of MHD LH2 void solution does not change 
very much with a magnetic parameter h > 0. The density pro- 
files also show shell-type morphology and the velocity still in- 
creases with radius (see Fig|3]l. Compared with non-magnetized 
cases shown in Figure|2]with the same void boundary x* and the 
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Figure 3: Global MHD LH2 void shock solutions with parameters (n = 
0.75, q = -0.5, y = 1.1875, h = 0.3) and a void boundary .v* = 1. The 
same format is adopted as in Figure [2] The void solution with K = 1 does 
not allow a magnetosonic shock and it merges into the MHD free-expansion 
solution. The void solution of K = 5 connects with solutions 1 and 2 whose 
parameters are: A = 12.325, B = 3.977, x si = 4, a si = 0.4915, v sd = 
2.2888, Xsa = 5.5577, a su = 0.1273, v 6u = 0.3533 (solution 1), and 
A = 18.390, B = 8.344, x sd = 7, a si = 0.2401, v sd = 4.2102, x m = 
8.6661, a su = 0.0941, v su = 3.2135 (solution 2), respectively. The void so- 
lution of K = 10 connects with solutions 3 and 4 whose parameters are: A = 
7.5520, B = 4.010, x sd = 4, a sd = 0.4478, v sd = 2.1273, x su = 4.2913, a su = 
0.1775, v sll = 0.8558 (solution 3), and A = 15.175, B = 7.706, x s(i = 7, a si = 
0.2317, v sd = 4.1652, x su = 8.0096, a m = 0.0999, v su = 3.1276 (solution 4), 
respectively. 



same parameter K, the peak density in the shell is lower and the 
shell width appears broadened in MHD cases. The K — 1 void 
solution in the non-magnetized case can involve a shock, oth- 
erwise it encounters the SCC. However, with the same value of 
K in the magnetized case, the void solution cannot harbor any 
shocks and merge into the MHD free-expansion solution defi- 
nitely. The K — 5 and K — 10 void solutions in the magnetized 
case can harbor MHD shocks. The MHD behaviour of the cor- 
responding upstream sides, from our numerical exploration, is 
all outflow (Solutions 1, 2, 3 and 4 of FigO. With a larger x S d, 
or a faster shock, the upstream outflow has larger velocity, both 
near the shock front and far from the centre. The shell-type 
appearance is commonly observed for LH2 void solutions, for 
hydrodynamic and MHD cases. 

We now consider MHD LH1 void solutions, for which the 
void boundary nx-v-0, a = is also a critical curve and the 
asymptotic solution approaching the void boundary is an eigen- 
solution. The regime of parameter in which the LH1 void solu- 
tion exists is q > 0, y+q > 1, h > (y > 1). Examples of MHD 
LH1 void solutions are shown in Fig [4] which do not encounter 
the critical curve and approach free-expansion asymptotic solu- 
tion (T37T > at large x. For this reason, the three velocity curves 
in the bottom panel of Fig |4] are nearly identical. The density 
and velocity increases as x increases. With a larger parameter q, 
the density profile appears more pronounced in the vicinity of 
void boundary. Because the free-expansion solution has infinite 
velocity and constant density far from the flow centre, the LH1 
void solutions would be more suitable for astrophysical model 
if they are matched with another branch of solution with finite 
velocity and density at large x by MHD shocks. Examples of 
shocks are shown in Fig [4] for q = 0.3 with downstream shock 
positions x s d = 3, x s( j = 5, and jt s d = 8, respectively. Again with 
a larger x s d, or a faster shock, the upstream outflow has a higher 
speed. 

Examples of MHD LH1 void solutions for the relativistic 
case of q — 2/3, y — 4/3 are displayed in Fig. [5] We set free pa- 
rameter C — 1 and the arbitrary parameter A on the shock to be 
A - I. With equations (O and (O, we obtain the MCC with 
a = 0.3897 and v = 0.015 Ix. In this case, the suitable expan- 
sion solution becomes the Einstein-de Sitter solution as shown 
in Fig. [5] Similarly, we are able to construct MHD shocks to 
match LH1 void solutions with another branch of solution ( 1321 
with finite velocity and density far from the flow centre. 

In general, the reduced density of MHD LH1 void solu- 
tions increases with increasing x, and the density near the up- 
stream side of shock is very low, as compared with the den- 
sity near the downstream side of shock (see Figs. [4] and E] for 
shocks). For the corresponding upstream solutions, the den- 
sity decreases as x increases and tends to zero at large x. With 
MHD shocks, we obtain again the shell-type morphology for 
the density. This shell-type morphology here is somehow dif- 
ferent from the shell-type morphology of LH2 void solutions 
(e.g. Figs. |2]and|3]l. The LH2 void solutions have density peaks 
near the void boundary by themselves: the density increases 
and then decreases with increasing x. However, the LH1 void 
solutions must involve MHD shocks to have shell-type density 
profiles and the peak density is located on the downstream side 
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Figure 5: Semi-complete global MHD LH1 void shock solutions with param- 
eters n = 0.75, q = 2/3, y = 4/3, h = 0.3, C = 1 and the void boundary 
x* = 1 and x* = 4. The same format as Figure [2] is adopted. The MCC is 
a = 0.3897, v = 0.0151jc. The inner void solution of X* = 1 connects with two 
upstream solutions whose parameters are: A = 0.2819, B = -0.3660, x sd = 
3, a si = 0.1577, v sd = 2.0388, x sll = 3, a sll = 0.0136, v su = -0.2542; 
and A = 2.0920, B = -0.7974, x si = 6, a si = 0.1666, v sd = 4.0166, x su = 
6, o; su = 0.0163, v 8u = -0.4443, respectively. The inner void solution of a:* = 4 
connects to one upstream solution whose parameters are: A = 2.7441, B = 
-1.5525, x sd = 8, a si = 0.1404, v sd = 5.5306, x su = 8, or 6u = 0.0097, v su = 
-0.7679. 



of shock. 



Figure 4: Semi-complete global MHD LH1 void solutions and shocks with the 
parameter (n = 0.75, h = 0.3) and the void boundary x* = 1. The same format 
as Figure|2]is adopted. The void solutions merge to the free-expansion solution 
)37t and the constant value of a to be a m = 0.238 (equation 1401 . The void 
solution of q = 0.3 connects with three upstream solutions whose parameters 
are: A = 0.0968, B = 1.6372, x sd = 3, a sd = 0.0299, v sd = 2.3340, x su = 
3.1374, a su = 0.0120, v su = 1.8694; A = 0.9024, B = 2.8925, ,v sd = 5, or sd = 
0.0774, v 8d = 3.7655, x su = 5.1942, a su = 0.0353, v 8U = 3.0033; and 
A = 5.0384, B = 5.2269, x sA = 8, a sd = 0.1300, v sd = 5.8366, x slI = 
8.2401, Q- su = 0.0689, v su = 4.7681, respectively. 



4. Hydrodynamic Self-Similar Void Solutions with a* ± 

By setting magnetic parameter h — 0, we readily obtain a 
group of self-similar nonlinear ODEs describing the hydrody- 
namics of a general polytropic gas with specific entropy con- 
served along streamlines. We shall choose a non-zero or* at 
the void boundary nx — V = 0, or parameter K, in constructing 
solutions with a considerable freedom. We also insert hydro- 
dynamic shocks to obtain semi-complete global solutions satis- 
fying the asymptotic condition that v(x — > oo) = 0. According 
to Table Q] asymptotic behaviours near the void boundary can 
be generally classified as Type-N (q = 0, n + 1), Type-N2 
(q = 0, n = 1) and Type-D (q > 0) separately. 

4.1. Cases of q = 

With q — and h = 0, the flow system is reduced to a con- 
ventional polytropic gas obeying n + y — 2 and have a constant 
specific entropy everywhere at all times. Such hydrodynamic 
flows are systematically and carefully analyzed and discussed 
by Wang & Lou (2007). Hu & Lou (2008a) constructed void so- 
lutions in a conventional polytropic flow to model the so-called 
'champagne flows' in H II regions. In Figure [6] we present 
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Figure 6: Void solutions with n = 0.85, y = 1.15, q = 0, h = for a con- 
ventional polytropic gas. The same format as Fig. [2] is adopted. The void 
solution of x* = I, a* =2 connects with solutions 1, 2 and 3 whose pa- 
rameters are: A = 2.2257, B = -1.8402, x sA = 2, a si = 1.4849, v sd = 
1.3202, x su = 2.5521, a su = 0.1870, v su = -1.6783 (solution 1); A = 
2.5682, B = 0, x si = 2.314, a si = 1.2452, v sd = 1.4646, x su = 2.4721, a slI = 
0.2797, v sll = -0.2876 (solution 2); and A = 6.4646, B = 3.2026, x si = 
3, Qf sd = 0.8228, Vsd = 1.7353, ,v su = 3.0056, a su = 0.5033, v 8U = 1.2203 
(solution 3), respectively. 

examples of void solutions with hydrodynamic shocks across 
the sonic singular surface. For the same void solution, by prop- 
erly choosing the downstream shock position x^, we can obtain 
various dynamic behaviours on the upstream side: outflow (e.g. 
B > in solution 3 of Figure |6j, inflow (e.g. B < in solu- 
tion 1 of Figure|6]l and contraction (e.g. B = in solution 2 of 
Figure |5J. Basically, upstream dynamic behaviours depend on 
the void boundary x*, the density at the void boundary a", and 
the downstream shock position x S( i- Extensive numerical explo- 
rations reveal that the isothermal case (Type-N2) is similar to 
other q — cases regarding the void solutions. Lou & Zhai 
(2009) provide a detailed analysis for isothermal voids. 

4.2. Cases of q > 

For q > 0, asymptotic solution behaviours at the void bound- 
ary are of Type-D. We can construct void solutions with shocks. 
Examples of such void shock solutions with different values of 
q are shown in Figs[7jand[8] 

Dynamic behaviours of void solutions depend on the void 
boundary x* and the density parameter K. Without encounter- 
ing the sonic critical point, some void solutions (see curves 2, 
3 and 4 in Fig [7J merge into one kind of asymptotic expan- 
sion solutions (e.g. asymptotic free-expansion solution, with 
V ~ 2x/3 for q < 2/3, Einstein-de Sitter solution for q = 2/3, 
or asymptotic thermal-expansion solution for q > 2/3). As 
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Figure 7: Inflow velocities -v(.v) of void solutions with n = 0.85, y = 1.2875, 
q = 0.5, h = 0. The upper dotted curve is the sonic critical curve (SCC), and 
the lower dotted curve is the void boundary nx — v = 0. The solid curve 1 
is a void solution with x* = 0.4, crossing the upper SCC smoothly at x = 4, 
and matching asymptotic solution (32) with A = 19.852 and B = 1.286. Solid 
curves 2, 3, 4 are integrated from x* = 1, x* = 2, x* = 3, with K = 0.29, 
and v' given by expression )16t . These solutions merge to the free expansion 
asymptotic solutions with form v = 2jc/3 + b, = 2/3. The void solution of 
x* = 1 and K = 0.29 connects with solutions 5 and 6 whose parameters are: A = 
3.241, B = -0.987, x sd = 2, a sd = 0.502, v sd = 1.246, x m = 4.045, a m = 
0.114, v su = -0.607 (solution 5), andA = 2.173, B = -0.143, jc sd = 2.4, a sd = 
0.435, = 1-479, ,v su = 3.425, a su = 0.126, v su = 0.136 (solution 6), 
respectively. 

shown by Figure [7] void solutions with different void bound- 
aries merge into the same asymptotic free-expansion solution, 
except for a slightly different b parameter. Again, we shall con- 
nect these solutions with another asymptotic solution of finite 
velocity and density at large x by shocks (see curves 5 and 6 
in Fig. |7J. By properly choosing the downstream shock posi- 
tion, we could readily obtain an inflow (v < 0) or an outflow 
(v > 0) for the upstream side of a shock. The radial flow ve- 
locity tends to zero at large radii. Some void solutions cross the 
critical curve smoothly (see curve 1 in Fig. |7J. 

As shown in Fig. [8] of a relativistically hot gas, the critical 
curve can be obtained analytically from equations (fT2l) and ( fT3l 
with a = 3.2843 and v = -2.1772x Void solutions merge 
into the Einstein-de Sitter solution as expected. The Einstein- 
de Sitter solution has a diverging velocity at large x. Again, this 
solution can be connected with an asymptotic solution (f32l > of 
finite velocity and density at large x via shocks. 

5. MHD Self-Similar Void Solutions with a* + 

The magnetic force may play a key role in the vicinity of void 
boundary and smooths out all divergence of the non-magnetized 
cases. In a MHD flow (q > 0), a* = at the void boundary is 
allowed for non-trivial void solutions with h > 0, e.g. MHD 
LH1 void solution. In this section, we compare a* — and 
a* > cases, and see if the boundary value of a* modifies 
characteristic behaviours of void solutions. 
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Figure 8: Void solutions and the Einstein-de Sitter solution with n = 0.85, 
y = 4/3, q = 2/3, h = and C = 1. Inflow velocities — v(x) are presented. The 
same format as FigurefJJis adopted. The critical curve is v = -2.1772a:. The 
solid curves are integrated from the void boundary at X* = 1, x* = 2, x* = 3, 
with K = 0.2 and v' given by expression 1 1 (it . which match with the Einstein-de 
Sitter solution at large x. The void solution with x* = 1, K = 0.2 connects to 
one upstream solution whose parameters are: A = 0.188, B = -12.131, JE s d = 
1.2, a si = 3.453, v sd = 0.874, x su = 1.2, ot su = 0.0524, v su = -8.566. 

5.7. Case of q = Ofor a conventional polytropic gas 

Examples of global MHD void solutions for the case of q — 
are shown in Fig. [9] Numerical computations show that along 
the MCC, velocity gradient V has one positive and another 
negative eigenvalues, corresponding to the MCC being sad- 
dle points (e.g. Jordan & Smith 1977). Solutions approach- 
ing this MCC, may either cross the critical curve smoothly and 
match with asymptotic solution (l32l > of finite velocity and den- 
sity at large x (see solution curve 1 in Fig. or be turned 
back smoothly to match with another branch of solutions and 
merge into the Einstein-de Sitter solution at large x (see solu- 
tion curves 2, 3 and 4 in Fig. By adjusting the a* value at 
the void boundary, we can make the solution crossing the MCC 
smoothly. The only difference to integrate curves 1 and 2 in Fig. 
[9] is the a* value. Compared with the case without magnetic 
field (see Fig. |6), a void solution in this case can merge into 
the Einstein-de Sitter solution at large x without encountering 
the MCC. Void solutions can be connected with outer branch of 
the asymptotic solutions of finite velocity and density by MHD 
shocks (see solution curves 5 and 6 in Fig. 01. Similarly by 
adjusting the downstream shock position x s d, one void solution 
on the downstream side can be connected to various upstream 
solutions with different behaviours at large x. 

5.2. Case of q > 

An example of q = 0.5 is shown in Fig. [10] which can be 
compared with Fig. [7] The case of q = 2/3 describes a rela- 
tivistically hot gas. MHD void solutions in such case is shown 
in Fig. [TTJfor a comparison with Fig. [8] 

Again the density at the void boundary a* can be set to ei- 
ther zero, which leads to the LH1 void solution with eigenso- 
lution dT3T > (see curves 2, 3 and 4 of Fig. [TO] and curve 1 of 
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Figure 9: MHD void solutions with n = 0.85, y = 1.15, q = and h = 0.3. 
Inflow velocities -v are plotted. The same format as Fig. 0is adopted. Solid 
curve 1 is a void solution with x* = 1 and or* = 5, crossing the MCC smoothly 
at x = 4 to match with asymptotic solution i32i with A = 32 and B = 5.413 at 
large x. Solid curves 2, 3, 4 are integrated from the void boundary at x* = 1, 
x* = 2, x* = 3 respectively all with a* = 2 and V given by expression \25\ . 
These solutions are limited by the upper MCC and merge into the Einstein-de 
Sitter solutions (v = 2jc/3 , a = const). The void solution of x* = 1 and 
a* = 2 connects with solutions 5 and 6 whose parameters are: A = 0.205, B = 
4.280, x sA = 3, a si = 0.556, v sd = 1.714, x m = 3.159, a su = 0.191, v su = 
0.123 (solution 5), and A = 0.169, B = 1.128, x si = 2.2, a s(i = 0.868, v sd = 
1.356, x su = 3.306, a su = 0.126, v su = -2.492 (solution 6), respectively. 



Fig. rm >. or a nonzero finite value, which leads to a Type-Nq 
behaviour of equations (f27t — <f30b (see curves 2', 3' and 4' of 
Fig. [10] and curves 2 and 2' of Fig. \TT[ . In the absence of a 
magnetic force, only the Type-D behaviour is allowed in this 
range of q and the density diverges on the void boundary. If not 
encountering the MCC, these solutions merge into one kind of 
expansion solutions at large x (i.e. free-expansion for q < 2/3, 
Einstein-de Sitter solution for q = 2/3, and thermal-expansion 
for q > 2/3). This property is the same as hydrodynamic cases, 
and we note that the constant density of the free-expansion 
solution does depend on the magnetic parameter h. This sug- 
gests that for such void solutions, the magnetic force plays an 
important role. Similar to non-magnetized cases, we can avoid 
the velocity divergence of the expansion solutions by matching 
such solutions with asymptotic solutions of finite velocity and 
density at large x via MHD shocks (e.g. curves 5, 6 and 7 of 
Fig. [10] and curves 3, 4 and 5 of Fig. ITTb . or making the void 
solutions crossing the MCC smoothly (see curve 1 of Fig. [10] 
and compare it with curve 1 of Fig. [7]i. From Figs. [10] and 
[TT1 the LH1 void solutions and the Type-Nq behaviour appear 
quite similar in terms of velocity profiles, except near the void 
boundary and different b parameter in the corresponding expan- 
sion solutions. We will discuss the influence of the initial mass 
density a* presently. 

The MCC has the form a = 3.0088, v = -2.2087x for the 
case of C = 1, y = 4/3, n = 0.85, q = 2/3 and h = 0.3. The be- 
haviour of MHD void solutions is similar to the hydrodynamic 



13 



n=0.85,q=2/3,h=0.3 



n=.85,q=0.5,h=0.3 



2 2 



-2 

>-A- 
-6 

-8- 
-10- 




012345678 

X 



9 10 




Figure 11: MHD void solutions with n = 0.85, y = 4/3, q = 2/3, h = 0.3 
and C = 1 . -v is plotted. The upper dotted curve is the MCC satisfying v = 
-2.2087a" and the lower dotted curve is the void boundary nx—v = 0. The solid 
curves 1 is LH 1 void solution, integrated from the void boundary at x* = 1 with 
a = 0. The solid curves 2 and 2' are Type-Nq void solutions, integrated from the 
void boundary at x* = 1 with a* = 2 and a* = 10. These solutions merge into 
the Einstein-de Sitter solution at large radii. The void solution 1 of x* = 1 , a* = 
connects with solutions 3 and 4 whose parameters are: A = 0.379, B = 
1.255, x sd = 3, a sd = 0.114, v sd = 2.095, x su = 3, a su = 0.0413, v su = 
1.292 (solution 3); and A = 0.848, B = 1.680, * sd = 4, a sd = 0.123, v sd = 
2.750, x su = 4, a su = 0.0459, v 8u = 1.661 (solution 4), respectively. The void 
solution 2' of x* = 1 and a* = 10 connects with solutions 5 whose parameters 
are: A = 2.699, B = -4.921, jc sd = 3, a si = 0.864, v sd = 1.544, x su = 
3, a su = 0.149, v su = -3.298. 



Figure 10: MHD void solutions with n = 0.85, y = 1.2875, q = 0.5, h = 0.3. 
Inflow velocities -v(x) are plotted. The same format as Fig. |7]is adopted. The 
solid curve 1 is a void solution with the void boundary at x* = 0.3, crossing 
the MCC smoothly at x = 5 to match with asymptotic solution (32) with A = 
25.977 and B = 1.147. The solid curves 2, 3, 4 are MHD LH1 void solutions, 
numerically integrated from the void boundary at x* = 1, x* = 2, x* = 3 
with a* = 0. Respectively, the dashed curves 2', 3', 4' are Type-Nq void 
solutions, integrated from the void boundary at x* = 1, X* = 2, x* = 3 with 
a* = 2. These solutions merge into the free expansion asymptotic solutions 
in the form of v = 2x/3 + b at large x with a m = 0.238. The void solution 
2 of x* = 1, of* =0 connects with solutions 5 and 6 whose parameters are: 
A = 0.256, B = 1.692, .v sd = 2, a sd = 0.0646, v sd = 1.475, x su = 3.662, a su = 
0.0182, v su = 1.648 (solution 5), and A = 0.805, B = 2.119, x si = 3, a s(l = 
0.0943, v sd = 2.142, .v su = 4.731, a su = 0.0299, v su = 1.991 (solution 6), 
respectively. The void solution 2' of x* = 1 and a* = 2 connects with solution 
7 whose parameters are: A = 2.610, B = 0.208, jr sd = 3, ci sd = 0.292, v sd = 
1.772, x su = 4.142, a su = 0.102, v su = 0.445. 



case. In such cases, it is unlikely for void solutions to encounter 
the MCC and the only way they cross this singular surface is via 
MHD shocks (see curves 3, 4 and 5 in Fig. fTTT i. 

We further investigate the influence of a* value at the void 
boundary on semi-complete global solutions. By comparing so- 
lutions from the same void boundary x* and different a* values 
(see curves 1, 2, 2' in Fig. [TTJ, it appears that with larger a* on 
the void boundary, the void solutions converge to the asymp- 
totic solution more slowly. This means that with larger density 
gradient on the void boundary, the system has a larger transi- 
tion zone where the magnetic force, the gravity and the thermal 
pressure force are all comparable. We refer to this zone as the 
void boundary layer. Outside the void boundary layer, the free 
expansion, thermal expansion or the Einstein-de Sitter solution 
would be a good approximation for the asymptotic behaviour. 
With an MHD shock inserted, the solution can be matched with 
an outflow (curves 3 and 4 in Fig. fTTb or an inflow (curve 5 in 
Fig. [TTJ. The dynamic behaviour of the outer upstream flow 
of the global void solution depends on the void boundary x*, 
the value of a* and the downstream shock point x s( j. From the 
same void boundary x*, we can adjust a* values to let the solu- 
tion either cross the MCC smoothly or merge into asymptotic 
expansion solutions (see curves 1 and 2 in Fig. |5). 

As a further discussion, we obtain void solutions with dif- 
ferent a* values at the void boundary with parameters (n = 
0.85, y = 1.7, q = 2, h = 0.3) (see Figure O- In such 
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Figure 12: Influence of the initial value of a* on the void solution. The 
five flow solution curves are integrated with n = 0.85, y = 1.7, q = 2, 
h = 0.3, from the void boundary at x* = 1, and different initial values of 
a* = 0.5, 0.4, 0.3, 0.2, 0.1 are chosen. The v{x) part of these five solutions 
are nearly identical (not shown here). From this figure, for x > 2.5, the five 
curves merge together, therefore the influence of the initial value of a* is only 
significant around 1 < x < 2.5. 

cases, there is no MCC. The solution curves are very similar 
and merge quickly as a single curve (i.e. asymptotic thermal 
expansion solution). This suggests that the general behaviour 
of the void solution is not influenced by the initial a* value. 
In other words, a* at the void boundary is fairly arbitrary and 
only influences the dynamical behaviour near the void bound- 
ary (e.g. void boundary layer). For a realistic astrophysical 
flow, the dynamics on the void boundary cannot be described in 
a self-similar manner, due to unavoidable diffusion processes. 
Therefore, choices of a* at the void boundary in our model 
serve only for starting a numerical integration. A more com- 
plete understanding of such system requires information for the 
relations between the a* value on the void boundary and the 
initial physical conditions that generate voids, such as density 
perturbations and growths in supernova explosions (e.g. Cao & 
Lou 2009) or hot fast winds in planetary nebulae (e.g. Lou & 
Zhai 2009). 

6. Astrophysical Applications 

Our self-similar solutions can be adapted to different astro- 
physical flow systems with various spatial and temporal scales. 
The sound parameter k determines the dimensional quantities 
in physical space and k varies for different astrophysical flow 
systems. From self-similar transformation ©, we have a rela- 
tion 

fcl-3<?/2 _ P 

(4nG)y- 1 G^(3n - 2)ipVM9 

= k JL (43) 

fi(4nGy- l Gi(3n-2)ipy- l Mi ' 

where T is the thermal temperature, kg is Boltzmann's constant, 
fi is the mean molecular (atomic) mass of gas particles and the 



second equality only holds for an ideal gas. Since the entropy is 
closely related to pp~ 7 , the increase of entropy across a shock 
front from the upstream side to the downstream side would lead 
an increase of k value in the same direction for q < 2/3 . For y > 
1 and q < 2/3, the temperature should also increase across a 
shock front in the same direction. It is still not trivial to estimate 
k values from relation d43l above, because the enclosed mass M 
varies in r. Numerical tests show that when q is not too large, 
setting q = does not influence the magnitude order of k value. 
Thus we use a simplified relation 



py(AnG)y- 1 npv- l (AnG)v- 1 ' 

This is identical with relation (59) of Lou & Wang (2006) for 
a conventional polytropic gas. The relation does depend on the 
value of y. For a late evolution phase of massive stars after 
the hydrogen burning, the central density and temperature are 
p c ~ 10 8 g cm" 3 and T c ~ 10 9 K. We estimate k ~ 10 16 - 10 17 
cgs unit, depending on the value of y. For the interstellar 
medium (ISM) in our own Galaxy, mainly composed of hydro- 
gen, p ISM ~ 10- 20 - 10- 26 g cm- 3 and r ISM ~ 10 - 10 6 K (e.g. 
Ferriere et al. 2001) and we estimate k ~ 10 9 - 10 24 cgs unit, 
depending on the value of y in the range ofl <y<4/3. 

The parameters we have adopted in our model are (n, y, q, h) 
with the relation y — 2 - n + (3n - 2)q/2. The physical meaning 
of these parameters is clear. Parameters (n, y, q) are relevant 
for general polytropic processes. By setting q — 0, we retrieve 
the conventional polytropic gas with a constant specific entropy 
everywhere at all times and require n + y — 2. The polytropic 
index y is an approximation commonly invoked when energetic 
processes are not known (e.g. Weber & Davis 1967). For ex- 
ample, when applying our model to an exploding stellar enve- 
lope, such as supernovae, y would be close to unity, indicating 
a tremendous energy deposit. To apply our model to a slowly- 
evolving ISM, y should be very close to ratio of specific heats 
Cp/c v for an adiabatic process. This is consistent with the dy- 
namic evolution shown by our self-similar solutions. For a fixed 
x value, the corresponding radius r expands with time obeying 
a power law of ~ ?". For n > 1, r expands faster and faster, 
implying a continuous energy input into a gas flow. Another 
role of n is that it scales the initial density distribution of a gas 
flow. According to asymptotic solution (1321 1. the mass density 
scales as x~ 2 ^ n at large x. The initial condition with t —> + cor- 
responds to the asymptotic boundary condition with x — > oo, 
so the scaling parameter n determines the initial density profile 
when the solution takes asymptotic form (l32l at large x. 

Our general shock void solutions may be adapted to model 
planetary nebulae. In the late stages of stellar evolution, the 
compact star becomes an intense source of hot fast stellar 
wind and photoionization. The fast wind catches up with the 
fully photoionized slow wind and supports a fast wind bub- 
ble of hot gas. Chevalier (1997) developed an isothermal self- 
similar model without self-gravity to study the expansion of a 
photoionized stellar wind around a planetary nebula (see also 
Meyer 1997). The key idea of Chevalier (1997) is that the in- 
ner edge of the slow wind forms a contact discontinuity with 
the stationary driving fast wind. We have shown that this con- 
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tact discontinuity in gravity-free cases corresponds to the void 
boundary in our formulation. Lou & Zhai (2009) presented an 
isothermal model planetary nebula involving an inner fast wind 
with a reverse shock; this shocked wind is connected to an ex- 
panding self-similar void solution through an outgoing contact 
discontinuity. In their model, the self-gravity is included and 
a variety of flow profiles are possible. We here provide a the- 
oretical model formulation in a more general framework with 
a polytropic equation of state and the inclusion of self-gravity 
and the magnetic force. Our MHD shocked void solutions are 
also suitable to describe the self-similar dynamics of planetary 
nebulae combined with effects of central stellar winds and pho- 
toionization. 

Another astrophysical context to apply our MHD void shock 
solutions is the expansion of H II regions surrounding new-born 
protostars, especially for "champagne flows" (e.g. Hu & Lou 
2008a). Ultraviolet photons from nascent nuclear-burning pro- 
tostars fully ionize and heat the surrounding gas medium and 
drive H II regions out of equilibrium. Such H II regions ex- 
pand and gradually evolve to a "champagne flow" phase with 
outgoing shocks. 

As a more detailed application of our void solutions, 
we revisit below the scenario for core-collapse supernovae, 
which has been investigated numerically over years (e.g. see 
Liebendorfer et al. 2005 for an overview). Neutrino-driven 
models are widely adopted for explaining the physical mech- 
anism of type-II supernovae. The core collapse and bounce 
create a tremendous neutrino flux and within several hundred 
milliseconds after the core bounce, neutrino sphere is largely 
trapped and deposit energy and momentum in the dense bary- 
onic matter. A typical scenario is that the neutrinos drive the 
stellar materials outwards, deposit large amount of outward mo- 
mentum and re-generate the delayed rebound shock to push out- 
wards. Janka & Muller (1995, 1996) successfully obtained nu- 
merical simulation for the first second of type-II supernovae 
based on such a scenario. Many recent numerical studies, with 
more careful consideration of the convection in the stellar enve- 
lope and diffusion processes, also confirm and consolidate the 
viability of such a neutrino reheating process (e.g. Buras et al. 
2006; Janka et al. 2007, 2008). From simulations of Janka & 
Muller (1995, 1996) on progenitor stars with a mass range of 
~ 8 - 15M , the neutrino sphere stops depositing energy ~ 0.5 
s after the core bounce and then decouples from the baryonic 
matter. Once decoupled, neutrinos quickly escape from the stel- 
lar interior and may leave a cavity between the centre and the 
envelope within the exploding progenitor star. At the centre of 
the cavity may lay a nascent neutron star, a stellar mass black 
hole, or even shredded debris (e.g. Cao & Lou 2009). 

We now show that the gravity of a remnant central object (if 
not completely destroyed by the rebound process) on the ex- 
panding stellar envelope may be neglected under certain situa- 
tions. The equivalent Bondi-Parker radius r^p is 



where M» is the mass of the central object and a is the sound 
speed of the surrounding medium, which mainly depends on 



temperature. Far beyond rep, the gravity of a central mass may 
be ignored. In the following, we will see that at ~ 1 s after the 
core bounce, the stellar envelope has a temperature of the or- 
der of ~ 10 8 K (see Fig. [T4l > and a corresponding sound speed 
squared a 2 ~ 10 17 cm 2 s~ 2 . With M, ~ M Q , the Bondi-Parker 
radius at 1 s is estimated to be of the order of ~ 10 8 cm which 
is roughly the same as the void boundary. Therefore at the be- 
ginning, the gravity of the central object is only marginally ig- 
norable. From self-similar transformation ©, the Bondi-Parker 
radius, which is proportional to 1 /a 2 ^ pip, has the time depen- 
dence of f 2 ~ 2 ". The void boundary expands as t". As long as we 
require n > 2/3, the Bondi-Parker radius expands slower than 
the void boundary; in other words, the gravity of the central ob- 
ject may be ignored shortly after the core bounce. We can then 
presume that the central cavity is an approximate void and ap- 
ply our self-similar void solutions. Similar approximation has 
been applied to other astrophysical flow systems, in which the 
gravity of central object can be neglected with respect to the 
dynamics of the surrounding gas (e.g. see Tsai & Hsu 1995; 
Shu et al. 2002; Bian & Lou 2005 and Hu & Lou 2008 for 
applications to shocked "champagne flows" in Hu regions). 

The model framework of self-similar dynamics described in 
this paper implies that the radius of a spherical shock r s evolves 
with time t in a power-law manner, i.e. r s = £ 1/2 x s f". We may 
regard k ' 2 x s and n as two free parameters and attempt to fit 
the self-similar evolution to the results of a numerical simula- 
tion by Janka & Muller (1996; e.g. their case 03c with relevant 
parameters specified in the caption of our Figure [T3l). We use 
k = 4 x 10 16 cgs unit for the inner part of the void solution (i.e. 
the downstream side of an outgoing shock). The best fit model 
is achieved at n = 1 .57 and the downstream shock position (or 
speed) of x s = 7.36 (see Fig. [T3l ). The self-similar evolution 
fits almost perfectly with the simulation. It is striking to obtain 
such a good agreement for shock evolution, because the numer- 
ical simulation of Janka & Muller (1996) employed an equation 
of state that contains contributions from neutrinos, free nucle- 
ons, a-particles and a representative heavy nucleus in nuclear 
statistical equilibrium. In other words, their model carries dis- 
tinct features. With our simple approximation, we essentially 
parameterized all these complicated energetic processes by a 
general polytropic equation of state. We expect to obtain dif- 
ferent best-fit scaling parameter n, in comparison with numer- 
ical results under different conditions, such as higher or lower 
initial neutrino luminosity. This fitting is suggestive that the 
polytropic approach is a fairly good approximation for shock 
evolution, and physically the rebound shock expands in a self- 
similar manner. 

The numerical simulation of Janka & Muller (1996) ends at 
~ 1 s. Within this duration, neutrinos deposit enough momen- 
tum and kinetic energy in a shocked stellar envelope and the 
star is set to explode as the rebound shock emerges from the 
stellar photosphere. The subsequent dynamic evolution, includ- 
ing the travel of the rebound shock, can be readily described by 
our self-similar model. Soon after neutrinos decouple from the 
baryon matter, no more energy is provided and indeed the sys- 
tem begins to lose energy by radiation processes. Thus, we can 
no longer apply n — 1.57 further. We therefore use their model 
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Figure 13: Shock positions as function of time after a core bounce. The solid 
curve shows the self-similar evolution, with parameters n = 1.57, k = 4 X 10 16 
cgs unit and the downstream shock position x s = 7.36. The asterisks show the 
result of numerical simulation by Janka & Miiller (1996) for a one-dimensional 
model of the core collapse of a progenitor with mass 15 M Q and the iron core 
mass 1.31 M , and the initial neutrino luminosity 2.225xl0 52 erg s~'. 

parameters at ~ 1 s as the initial input parameters of subsequent 
dynamic evolution. We found that if « = 0.8, y — 1.2, q = 
(i.e. a conventional polytropic gas), our self-similar model 
gives appropriate solutions and we do not include a random 
magnetic field in this preliminary illustration. The void solu- 
tions with these parameter are of Type-N and relevant exam- 
ples are also shown in Fig. [6] The cavity radius is taken to be 
~ 1000 km at t — Is and the corresponding void boundary is 
then x* = 0.5. From the simulation case 03c of Janka & Miiller 
(1996), the rebound shock is at ~ 1.3 x 10 9 cm at t = 1 s, which 
correspond to a downstream shock position x s a = 6.5. We note 
that for these parameters the mass density cannot be set to zero 
at the void boundary, hence we should choose a* properly such 
that the mass density on the void boundary at t — 1 s is equal 
to the mass density given by the case 03c. The void solution at 
t = 1 s is shown in Fig. [14] 

The solution shown in Figure [l4]corresponds well to the en- 
velope of an exploding massive star. The enclosed mass is 
~ 25M with a radius of ~ 10 12 cm, grossly consistent with 
typical masses and radii of O and B stars. The enclosed mass 
mainly depends on the density a* on the void boundary; by 
regarding a* as a free parameter, the self-similar dynamics is 
capable of modelling stars with different masses. The solution 
shows an expansion velocity ~ 10 9 cm s _1 , a typical expan- 
sion velocity for type-II supernovae. The temperature increases 
from the void boundary to the rebound shock and decreases 
with r outside the rebound shock. The temperature rises to 
~ 10 8 K, grossly consistent with typical supernova temperature. 
Furthermore, we can compute the total energy of the entire sys- 
tem. We should consider the kinetic energy E^, gravitational 
energy E g and the internal energy Ej. When calculating the in- 
ternal energy, we assume that the gas is single-atomic, whose 
degree of freedom is 3. At t — Is, the total energy E tota i is 




[cm] 



Figure 14: Self-similar void solution for an exploding progenitor star at t = Is. 
From top to bottom, the panels show the density, radial flow velocity, thermal 
pressure, enclosed mass and the thermal temperature. The solution is obtained 
with parameters n = 0.8, y = 1.2, q = 0, h = 0, k = 4 X 10 16 cgs unit, 
void boundary x* = 0.5, density at the void boundary a* = 0.001 and the 
downstream shock position x s d = 6.5. Correspondingly at the moment shown, 
the void boundary is at 10 s cm and the rebound shock is at 1.3 X 10 9 cm. Our 
solution is numerically integrated until 10 12 cm. 
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6.54 x 10 50 erg, in which the kinetic energy E k - 1.63 x 10 s 1 
erg, the gravitational energy E g = -1.10 x 10 51 erg, and the 
internal energy £, = 1.21 x 10 50 erg. We see that the kinetic 
energy is the major energy source, and the radial motion of 
the fluid has not been dissipated much to random motions, in 
agreement with the exploding star scenario. The total energy of 
our self-similar model is well consistent with the total energy 
9.5 x 10 50 erg given by the Case 03c of simulation (Janka & 
Miiller 1996) which suggests that the energy of the stellar enve- 
lope of the self-similar model comes from the neutrino process 
before 1 s. Our calculations reveal that the energy does not vary 
much with time, consistent with our assumption that the system 
is nearly adiabatic. We conclude that in general, our self-similar 
void solution is plausible as it shows typical value of velocity, 
pressure, temperature, enclosed mass and total energy. We can 
also adjust the self-similar parameters to produce solutions for 
different objects. 

Here we discuss two immediate utilizations of our self- 
similar solutions. First, as long as we know the evolution of a 
rebound shock, we can estimate the time when the shock breaks 
out of the stellar envelope. Assuming the stellar photosphere at 
a radius of ~ 10 12 cm, with the relation r s h oc k K t", we estimate 
that the shock travels to the photosphere at ~ 4 x 10 3 s (i.e. 
> 1 hour) after the core rebounce. From then on, we should 
be able to detect the massive star in act of an explosion. With 
advanced instrument, astronomers are now detecting more and 
more shock breakout events, and consolidate the association of 
long y-ray bursts (GRBs) and supernovae (e.g. Campana et al. 
2006). From Figure [14] the temperature around a shock in the 
stellar envelope is in the range of ~ 10 s K and gradually de- 
creases to the order of ~ 10 7 K as the shock breaks out of the 
stellar atmosphere. Such a temperature range will give rise to 
X-ray radiations. With our dynamical rebound shock model, 
coupled with radiation (e.g. the thermal bremsstrahlung) and 
transfer processes, we can calculate early X-ray emissions from 
supernovae in act of an explosion (Lou & Zhai 2009), and in 
turn, we may infer properties of GRB/SN progenitor by ob- 
servations of shock breakout diagnostics. Hu & Lou (2008b) 
presented some preliminary model calculations along this line 
and found sensible agreement with X-ray observations of SN 
2008D (e.g. Soderberg et al. 2008; Mazzali et al. 2008). 

Secondly, after long temporal lapses, we intend to relate the 
central cavity of the stellar envelope with hot bubbles observed 
in our own Galaxy. A typical SNR grows for ~ 1 .5 Myr and 
reaches a radius of ~ 50 pc (e.g. Ferriere 2001). We assume 
that the cavity and the stellar envelope (gradually evolving into 
a SNR) expand in a self-similar manner as the solution shown in 
Figure [14] after a SN breakout event. The cavity radius expands 
as f, and 1.5 Myr later, the radius of the cavity becomes ~ 3 
pc. Possible explanations for the discrepancy are: first, in the 
ISM the sound scaling factor k is much larger than that in the 
stellar envelope; therefore we can no longer assume k to be an 
overall constant. Secondly, the typical scale of a SNR includes 
both the cavity radius and the radius of the matter shell in the 
surrounding; and the gas shell can well spread out to the ISM 
since the pressure in the ISM is very low. 

As the last part of our discussion, we emphasize the shell- 



type solutions such as LH1 void solutions and LH2 void so- 
lutions, and corresponding examples are shown in Figs. [2] [3] 
[4] and [5] These void solutions have zero density at the void 
boundary and thus a sensible continuity across the void bound- 
ary. It is essential to consider the magnetic field for the LH1 
void solutions, as our study shows that the magnetic field is 
indispensable to obtain shell-type LH1 solutions. As the shell- 
type SNRs are commonly observed, and the magnetic field gen- 
erally exists in SNRs, our self-similar void solutions are likely 
to be a sensible dynamic approximation. In cases with q > 2/3, 
LH1 shell-type solutions merge into the asymptotic thermal- 
expansion solution at large x with a divergent velocity. These 
solutions do not encounter the MCC, indicating that the fluid 
keeps sub-magnetosonic for the entire flow system and such 
shell-type solutions cannot be matched with another asymptotic 
solutions of finite velocity and density at large x by an MHD 
shock. In reality, the outer layer of shell-type SNRs is bounded 
by the ISM. 

We define a shell width as the distance from the void bound- 
ary to the place where the density is e~ l of the peak density and 
perform numerical exploration to examine how the shell width 
depends on magnetic field strength h and the self-similar pa- 
rameter q. One example of LH2 solutions with n = 0.75, x* 
I , K — 1 is shown in Figure [15] For certain q, the shell width 
has a maximum value with h < 1 . The shell width increases 
with h rapidly in the regime h <§c 1 and decreases with h gradu- 
ally in the regime of h > 1 . For a larger q value, the maximum 
shell width is larger, and the shell width at large h is smaller. 
In general, we find that the magnetic field and the entropy dis- 
tribution (or parameter q) influence the shell-type morphology 
significantly. This implies that shell structures in supernova 
remnants may reveal some information on the magnetism of 
interstellar media or the internal energy input or output of the 
gas. 

7. Summary and Conclusions 

In a general polytropic MHD formulation, we obtain novel 
self-similar void solutions for a magnetofluid under self-gravity 
and with quasi-spherical symmetry. For our MHD void solu- 
tions, the enclosed mass within the void boundary is zero and 
our solutions are valid from the void boundary outwards. We 
have carefully examined MHD and hydrodynamic behaviours 
near the void boundary and obtained complete asymptotic solu- 
tions in various parameter regimes. For clarity, we consider the 
situation of density being zero at the void boundary and other 
situations separately. 

For the case of a* = at the void boundary, we derive novel 
LH1 solutions with q > 0, for which the void boundary is also 
a critical curve, and LH2 solutions with q < 0, for which the 
thermal pressure force becomes dominant approaching the void 
boundary. A random magnetic field must be present in order 
to construct LH1 solutions. Both LH1 and LH2 solutions have 
shell-type morphology in the density profile, whose peak den- 
sity and the shell width are mainly determined by the values of 
the magnetic parameter h and the parameter q. The shell width 



18 



n=0.75, x*=1, K=1 (LH2 solutions) 

5.5 i 1 1 1 1 1 1 1 1 




01 23456789 10 

h 

Figure 15: Shell width variation of LH2 solutions versus the magnetic param- 
eter h for different values of q parameter. The shell width is defined as the 
distance from the void boundary to the place where the density is e _I of the 
peak density. The common parameters for these solutions are n = 0.75, x* = 
1, K= 1. 

of the density profile expands with time also in a self-similar 
manner. 

For the case of a* at the void boundary, the situation is 
quite different between the hydrodynamic and MHD cases. In 
hydrodynamic cases, asymptotic behaviours in the vicinity of 
void boundary depends on parameter q and can be classified as 
Type-N (q = 0) and Type-D (q > 0). For the Type-D behaviour, 
the density at the void boundary diverges, and local diffusion 
process should occur. We systematically examined all these 
possibilities and present a few numerical solution examples. 
Our solutions are well compatible with previous self-similar so- 
lutions; for example, by setting q = 0, our formulation reduces 
to self-similar solutions for a conventional polytropic gas (Lou 
& Wang 2007). Without magnetic field, for both Type-N and 
Type-D behaviours, the thermal pressure force dominates at the 
void boundary; actually Type-N solutions can be regarded as a 
natural extension of LH2 solutions in the regime of q = 0, and 
Type-D solutions as a natural extension of LH2 solutions in the 
regime of q > 0. Therefore, for all sensible void solutions in 
hydrodynamic framework, the thermal pressure is the dominant 
force at the void boundary. By including a random magnetic 
field, which is ubiquitous in astrophysical plasmas, the diver- 
gence at the void boundary appearing in the type-D behaviour 
can be removed, and the dominant force on the void boundary 
becomes the magnetic force. 

Void solutions may go across the critical curve either 
smoothly or by an MHD shock and then merge into asymp- 
totic solution ( 1321 of finite density and velocity at large x. If 
void solutions do not encounter the critical curve, they gener- 
ally merge into one kind of asymptotic expansion solutions at 
large x with the velocity proportional to the radius. For q < 2/3, 



void solutions merge into asymptotic free-expansion solutions, 
for which the thermal pressure force is negligible. For q > 2/3, 
void solutions merge into asymptotic thermal-expansion solu- 
tions, for which the thermal pressure force is dominant. For 
q = 2/3, void solutions merge into the Einstein-de Sitter so- 
lution, a semi-complete global exact solution. We are free to 
choose the position of void boundary x*, the density on the 
void boundary a* (or the density parameter K) and the down- 
stream shock position x s d to construct various void solutions 
with different asymptotic dynamic behaviours far from the void 
centre, including inflows, outflows, contraction and breeze for 
upstream solutions. 

In this paper, we briefly discussed the case of q = 2/3 and 
thus y — 4/3 for a relativistically hot gas (Goldreich & Weber 
1980; Lou & Cao 2008; Cao & Lou 2009). One more parameter 
appears in the self-similar form of the equation of state denoted 
as C in this case. We show solution examples of C = 1. A more 
detailed study, including cases of C + 1 is forthcoming. 

Finally, we provide examples of applications of our self- 
similar MHD void solutions. In principle, our solutions can be 
adapted to various astrophysical plasmas with a central cavity. 
The scale of a system can vary within the upper limit of neglect- 
ing the universe expansion. The input and output of energy can 
be approximated by properly choosing parameters (n, y, k). 
As more self-consistent solutions of an MHD problem usu- 
ally require tremendous computational effort, our self-similar 
approach is valuable in conceptual modelling and in checking 
simulation results. We provide an application of our void so- 
lutions to the neutrino reheating mechanism for core-collapse 
supernovae and compare the dynamical results of self-similar 
solutions with previous numerical simulations. We find that our 
simplified model fits well with numerical simulations, suggest- 
ing that the self-similar approach is plausible and the explod- 
ing stellar envelope and the rebound shock do evolve in a self- 
similar manner. More specifically, we estimate that a rebound 
shock breaks out from the stellar photosphere ~ 4000 s after 
the core bounce, for a progenitor star of mass 25 M and radius 
~ 10 12 cm. We expect that, our general polytropic self-similar 
MHD solutions, coupled with radiative transfer processes, may 
offer physical insight for the rebound shock evolution of mas- 
sive stars as well as on GRB-supernova associations. We also 
indicate potential applications of LH1 and LH2 void solutions, 
with shell-type morphology, on the shell-type supernova rem- 
nants as well as hot bubbles in the interstellar medium. 
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A. MHD shock jump conditions 

MHD shocks can be constructed for self-similar solutions to 
cross the MSS. Despite discontinuities in pressure, mass den- 
sity, temperature, magnetic field and velocity across the shock 
front, we require conservations of mass, radial momentum, 
MHD energy across a shock front in the comoving shock refer- 
ence framework, respectively. They are 



[p(u s - u)]i = , 



p + p{u s -u) + 



<Bf> 



8tt 



= 0, 



p(u s - u) yp(u s — u) < B t > 
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(u s — u) 
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where u s is the shock travel speed in the laboratory frame- 
work of reference. Since we only consider the dominant trans- 
verse magnetic field parallel to the shock front in our theoretical 
model framework, the magnetic induction equation can be writ- 
ten as 



[(u s - uf < B] > ]j = . 



(49) 



The magnetic field average is made over a layer between r and 
r + dr, and such average can still describe the discontinuity in 
the radial direction. Strictly speaking, the magnetic fields have 
weak radial components normal to the shock front. We presume 
that such radial components are extremely weak compared to 
the transverse components on large scales. In common with the 
conventional shock analysis (e.g. Zel'dovich & Raizer 1966, 
1967), we use a pair of square brackets outside each expres- 
sion enclosed to denote the difference between the upstream 
(subscript '1') and downstream (superscript '2') sides across a 
MHD shock front. Note that the definitions of the downstream 
and upstream sides are in the reference framework where the 
shock front is at rest, and the specific entropy increases from 
the upstream side to the downstream side. The sound parameter 
k in transformation © is related to the polytropic sound speed 
and changes across a shock. We therefore relate upstream k\ 
and downstream &2 with a ratio factor A such that 



A 2 ki 



h\ = hi 



x\ — Axi 



(50) 
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The latter two in condition ( f50b are for the jump of the mag- 
netic energy density < B 2 > and the continuity of shock radius 
r s across the MHD shock front, which means while the dimen- 
sionless shock front x s has different values across the shock, 
they correspond to the same shock front radius r s . With the 
definition of h, h\ = h 2 actually means the transverse mag- 
netic field is proportional to the mass density across the shock 
front, consistent with Dyson & Williams (1997). Using h in- 
stead of h\ and h 2 for the magnetic parameter and substituting 
equation d50b into MHD shock conditions (|46]|-(|491, we obtain 
self-similar MHD shock conditions 

a\(nx\ - Vi) = Aa 2 (nx 2 - v 2 ) , (51) 



Ca 2-n + 3n q /2 x 2 q(nxi _ ^ 

, ha\x\ 
+a\(nx\ — Vi) + 



A 2 



2 

C a 2 -" +3 " q/2 x 2 2 q (nx 2 - vif 



+a 2 {nx 2 - v 2 ) + 



2 

2 , 2 r 

si l-n+3nq/2 

(r-D 



(52) 



(nx\ — vi) + ~Ca 1 



xxj 9 (njci — v\) q + 2ha\x\ 



= A 2 



/ \2 , _ l- n+ 3no/2 

(nx 2 - v 2 ) + —Cor, 



xx 2q (nx 2 — v 2 ) q + 2ha 2 x\ 



(r-D 

(53) 



in terms of the dimensionless reduced variables. Once we have 
(x 2 , a 2 , v 2 ) on the downstream side of a shock, we can deter- 
mine (jci, a\, vi) on the upstream side using MHD shock con- 
ditions (IBTl - (|53T > (Wang & Lou 2008) or vice versa. In cases 
of q — 2/3, there are only two independent relations among 
d5Tl>-(l53l (Lou & Cao 2008); therefore one may choose pa- 
rameter A > fairly arbitrarily, e.g. A = 1. 
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